/*
 * block.h
 *
 * Copyright (c) 2009, Eugen Stoian <stoian.e@gmail.com>
 * 
 * This library is free software; you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the Free
 * Software Foundation; either version 2.1 of the License, or (at your option)
 * any later version.
 * 
 * This library is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
 * details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
 * http://www.fsf.org/licensing/licenses/lgpl.txt
 */

#ifndef BLOCK_H_
#define BLOCK_H_
#include <vector.h>
namespace pxx
{
//!Class template used as base class for matrix class template
template<typename T>
class matrix
{
protected:
	void allocate();
	void initialize(int rows, int columns);
	void free();
public:
	matrix(int rows = 0, int columns = 0);
	matrix(matrix const &);
	matrix& operator=(matrix const &);
#ifdef __GXX_EXPERIMENTAL_CXX0X__
	matrix(matrix &&);
	matrix& operator=(matrix &&);
#endif

	virtual ~matrix();

	T operator()(int i, int j) const;

	T& operator()(int i, int j);
	/*!
	 * \fn operator+(matrix const & m) const
	 * \brief Addition operator
	 * \param m a const reference to a matrix object of the same size
	 * \throw 0 is thrown if matrices have incompatible sizes
	 * \sa matrix::add()
	 */
	inline matrix operator+(matrix const & m) const {
		return add(m);
	}
	/*!
	 * \fn operator-(matrix const & m) const
	 * \brief Substraction operator
	 * \param m a const reference to a matrix object of the same size
	 * \throw 0 is thrown if matrices have incompatible sizes
	 * \sa matrix::sub()
	 */
	inline matrix operator-(matrix const & m) const {
		return sub(m);
	}
	/*!
	 * \fn operator*(matrix const & m) const
	 * \brief Multiplication operator
	 * \param m a const reference to a matrix object with the vector count equal with the current object column count
	 * \throw 0 is thrown if matrices have incompatible sizes
	 * \sa matrix::mul()
	 */
	inline matrix operator*(matrix const & m) const {
		return mul(m);
	}
	/*!
	 * \fn operator*(vector<T> const & v) const
	 * \brief Multiplication operator
	 * \param m a const reference to a matrixps object with the vector count equal with the current object column count
	 * \throw 0 is thrown if matrices have incompatible sizes
	 * \sa matrix::mulv()
	 */
	inline vector<T> operator*(vector<T> const & v) const {
		return mulv(v);
	}
	/* \fn operator*(T d) const
	 * \brief Multiplication  to the right with a scalar operator
	 * \param d a T precision float
	 * \sa matrix::muls()
	 */
	inline matrix operator*(T d) const {
		return muls(d);
	}
	/*!
	 * \fn operator~() const
	 * \brief Transpose operator.
	 * \sa matrix::transpose()
	 */

	inline matrix operator~() const {
		return transpose();
	}

	/*!
	 * \fn rowCount() const
	 * \brief Returns the rows number of the matrix.
	 */
	inline int rowCount() const {
		return _rowCount;
	}
	/*!
	 * \fn columnCount() const
	 * \brief Returns the columns number of the matrix.
	 */
	inline int columnCount() const {
		return _columnCount;
	}

	matrix add(matrix const & m) const;

	matrix sub(matrix const & m) const;

	matrix mul(matrix const & m) const;

	matrix muls(T d) const;

	vector<T> mulv(vector<T> const & v) const;

	matrix transpose() const;

	matrix getBlock(int i0, int i1, int j0, int j1) const;

	vector<T> getRow(int r, int i0, int j0) const;

	vector<T> getColumn(int c, int i0, int j0) const;

	inline bool isSized() const{
		return _rowCount*_columnCount;
	}
	inline bool isAllocated() const{
			return  _data;
	}
protected:
  T ** _data;
	int _rowCount;
	int _columnCount;
};


template<typename T>
void matrix<T>::allocate(){
	_data = (T**) pxx::memalloc(_rowCount * sizeof(T*));
	for (int i = 0; i < _rowCount; ++i)
		_data[i] = (T*) pxx::memalloc(_columnCount * sizeof(T));
}

template<typename T>
void matrix<T>::initialize(int rows, int columns){
	_rowCount = rows;
	_columnCount = columns;
	if(isSized())
		allocate();
	else{
		_rowCount = 0;
		_columnCount = 0;
	}
}

template<typename T>
void matrix<T>::free(){
	if(_data) {
		for (int i = 0; i < _rowCount; ++i)
			pxx::memfree(_data[i]);
		pxx::memfree(_data);
	}
}

/*!
 * \fn matrix<T>::matrix(int rows, int columns)
 * \brief Constructor
 * \param rows positive integer
 * \param columns  positive integer
 * Constructs a matrix object of size rows x columns.
 * If rows and/or columns are 0 (or called with default arguments)
 * matrix object is uninitialized; in this case initialization is performed
 * at first assignment to an matrix object which is initialized.
 */
template<typename T>
matrix<T>::matrix(int rows, int columns) :
	_data(0), _rowCount(rows), _columnCount(columns) {
	std::cout<<"Constructor called"<<std::endl;
	if (isSized()) {
		allocate();
		for(int i =0; i < _rowCount; ++i)
			for (int j = 0; j < _columnCount; ++j)
				_data[i][j] = 0;
	}
	else {
		_rowCount = 0;
		_columnCount = 0;
	}
}

/*!
 * \fn matrix<T>::matrix(matrix const & m)
 * \param m a const reference to a matrix object
 * \brief Copy constructor
 */
template<typename T>
matrix<T>::matrix(matrix const & m) :
	_data(0), _rowCount(m._rowCount), _columnCount(m._columnCount) {
	std::cout<<"Copy constructor called"<<std::endl;
	if(isSized()){
		allocate();
		for(int i = 0; i < _rowCount; ++i)
			memcpy(_data[i],m._data[i],_columnCount*sizeof(T));
	}
	else {
		_rowCount = 0;
		_columnCount = 0;
	}
}

/*!
 * \fn matrix<T>::operator=(matrix const & m)
 * \param m a const reference to a matrix object of the same size. If current object is unsized m can have any size
 * \throw 0 is thrown if  m have different size.If current object is uninitialized m can have any size.
 * \brief Assignment operator
 * \throw 0 if current object and m have different sizes (rows x columns)
 */
template<typename T>
matrix<T>& matrix<T>::operator=(matrix const & m) {
	std::cout<<"Copy assignment called"<<std::endl;
	if(this!=&m){
		if(isAllocated()){
			if((_rowCount!=m._rowCount)||(_columnCount!=m._columnCount))
				throw 0;
				for(int i = 0; i < _rowCount; ++i)
						memcpy(_data[i],m._data[i],_columnCount*sizeof(T));
		}
		else {
			if(m.isAllocated()){
				initialize(m._rowCount, m._columnCount);
				for(int i = 0; i < _rowCount; ++i)
					memcpy(_data[i],m._data[i],_columnCount*sizeof(T));
			}
		}
	}
	return *this;
}

#ifdef __GXX_EXPERIMENTAL_CXX0X__
/*!
 * \fn matrix<T>::matrix(matrix && m)
 * \param m a right reference to a matrix object of the same size.
 * \brief Move constructor
 * \throw 0 if current object and m have different sizes (rows x columns)
 */
template<typename T>
matrix<T>::matrix(matrix && m)
:_data(m._data),_rowCount(m._rowCount),_columnCount(m._columnCount){
	std::cout<<"Move constructor called"<<std::endl;
	m._data = 0;
}
/*!
 * \fn matrix<T>::operator=(matrix && m)
 * \param m a right reference to a matrix object of the same size. If current object is unsized  m can have any size
 * \throw 0 is thrown if m have different size.If current object is uninitialized  m can have any size.
 * \brief Move assignment operator
 */
template<typename T>
matrix<T>& matrix<T>::operator=(matrix && m) {
	std::cout<<"Move assignment called"<<std::endl;
	if(isAllocated()){
			if((_rowCount!=m._rowCount)||(_columnCount!=m._columnCount))
				throw 0;
			free();
	}
	_rowCount = m._rowCount;
	_columnCount = m._columnCount;
	_data = m._data;
	m._data = 0;
	return *this;
}
#endif

/*!
 * \fn matrix<T>::~matrix()
 * \brief Destructor.
 */
template<typename T>
matrix<T>::~matrix() {
	free();
}

/*!
 * \fn matrix<T>::operator()(int i,int j) const
 * \param i integer greater than 0
 * \param j integer greater than 0
 * \brief Returns the value of the element at vector \m i and column \m j in the matrix;
 */
template<typename T>
T matrix<T>::operator()(int i, int j) const {
	return _data[i][j];
}

/*!
 * \fn matrix<T>::operator()(int i,int j)
 * \param i integer greater than 0
 * \param j integer greater than 0
 * \brief Returns a reference to the element at vector \m i and column \m j in the matrix;
 */
template<typename T>
T& matrix<T>::operator()(int i, int j) {
	return _data[i][j];
}

/*!
 * \fn matrix<T>::add(matrix  const & m) const
 * \brief Addition operation
 * \param m a const reference to a matrix object of the same size
 * \throw 0 is thrown if matrices have incompatible sizes.
 * Returns the sum of the current matrix with the matrix \m m
 */
template<typename T>
matrix<T> matrix<T>::add(matrix const & m) const {
	if ((_rowCount != m._rowCount) || (_columnCount != m._columnCount))
		throw 0;
	matrix result(_rowCount, _columnCount);
#ifdef _OPENMP
	if(omp_get_thread_num()!=0) { //Avoid nesting threads
#endif
	for (int i = 0; i < _rowCount; ++i)
		core::VAddV(_columnCount, _data[i], m._data[i], result._data[i]);
#ifdef _OPENMP
}
else {
#pragma omp parallel for
	for(int i = 0; i < _rowCount; ++i)
	core::VAddV(_columnCount, _data[i], m._data[i], result._data[i]);
}
#endif
	return result;
}

/*!
 * \fn matrix<T>::sub(matrix  const & m) const
 * \brief Substraction operation
 * \param m a const reference to a matrix object of the same size
 * \throw 0 is thrown if matrices have incompatible sizes.
 * Returns the difference between the current matrix and the matrix \m m
 */
template<typename T>
matrix<T> matrix<T>::sub(matrix const & m) const {
	if ((_rowCount != m._rowCount) || (_columnCount != m._columnCount))
		throw 0;
	matrix result(_rowCount, _columnCount);
#ifdef _OPENMP
	if(omp_get_thread_num()!=0) { //Avoid nesting threads
#endif
	for (int i = 0; i < _rowCount; ++i)
		core::VSubV(_columnCount, _data[i], m._data[i], result._data[i]);
#ifdef _OPENMP
}
else {
#pragma omp parallel for
	for(int i = 0; i < _rowCount; ++i)
	core::VSubV(_columnCount, _data[i], m._data[i], result._data[i]);
}
#endif
	return result;
}

/*!
 * \fn matrix<T>::mul(matrix const & m) const
 * \brief Multiplication operation
 * \param m a const reference to a matrixps object with the vector count equal with the current object column count
 * \throw 0 is thrown if matrices have incompatible sizes
 * Returns the multiplication of the current matrix with the matrix \m m
 */
template<typename T>
matrix<T> matrix<T>::mul(matrix const & m) const {
	if (_columnCount != m._rowCount)
		throw 0;
	matrix result(_rowCount, m._columnCount);
#ifdef _OPENMP
	if(omp_get_thread_num()!=0) { //Avoid nesting threads
#endif
	for (int j = 0; j < _columnCount; ++j)
		for (int i = 0; i < _rowCount; ++i)
			core::VMulSAddV(m._columnCount, m._data[j], _data[i][j], result._data[i], result._data[i]);
#ifdef _OPENMP
	}
	else {
#pragma omp parallel for
		for(int j = 0; j < _columnCount; ++j)
			for(int i = 0; i < _rowCount; ++i)
				core::VMulSAddV(m._columnCount,m._data[j],_data[i][j],result._data[i],result._data[i]);
	}
#endif
	return result;
}

/*!
 * \fn matrix<T>::muls(T d) const
 * \param d a T precision float
 * \brief Multiplication to the right with a scalar operation
 * Returns the multiplication of the current matrix with the scalar \m d
 */
template<typename T>
matrix<T> matrix<T>::muls(T d) const {
	matrix result(_rowCount, _columnCount);
#ifdef _OPENMP
	if(omp_get_thread_num()!=0) { //Avoid nesting threads
#endif
	for (int i = 0; i < _rowCount; ++i)
		core::VMulS(_columnCount, _data[i], d, result._data[i]);
#ifdef _OPENMP
	}
	else {
#pragma omp parallel for
		for(int i = 0; i < _rowCount; ++i)
			core::VMulS(_columnCount, _data[i], d, result._data[i]);
	}
#endif
	return result;
}

/*!
 * \fn matrix<T>::mulv(vector<T> const & v) const
 * \brief Multiplication to the right  with a vector
 * \param v a const reference to a vector object with the size equal with the current object column count
 * \throw 0 is thrown if objects have incompatible sizes
 * Returns the multiplication of the current matrix with the vector \m v
 */
template<typename T>
vector<T> matrix<T>::mulv(vector<T> const & v) const {
	vector<T> result(_columnCount);
	for (int i = 0; i < _rowCount; ++i) {
		result(i) = core::VDotV(_columnCount, _data[i], v._data);
	}
	return result;
}

/*!
 * \fn matrix<T>::transpose() const
 * \brief Transpose operation
 * Returns the transposed matrix.
 */
template<typename T>
matrix<T> matrix<T>::transpose() const {
	matrix result(_columnCount, _rowCount);
#ifdef _OPENMP
	if(omp_get_thread_num()!=0) { //Avoid nesting threads
#endif
	for (int i = 0; i < _rowCount; ++i)
		for (int j = 0; j < _columnCount; ++j)
			result._data[j][i] = _data[i][j];
#ifdef _OPENMP
	}
	else {
#pragma omp parallel for
		for (int i = 0; i < _rowCount; ++i)
			for (int j = 0; j < _columnCount; ++j)
				result._data[j][i] = _data[i][j];
}
#endif
	return result;
}

/*!
 * \fn matrix<T>::getBlock(int i0, int i1, int j0, int j1) const
 * \param i0 a positive integer not greater than matrix<T>::rowCount()-1
 * \param i1 a positive integer not greater than matrix<T>::rowCount()-1
 * \param j0 a positive integer not greater than matrix<T>::columnCount()-1
 * \param j1 a positive integer not greater than matrix<T>::columnCount()-1
 * \return a matrix<T> object
 * \brief Returns the copy of the block starting at (i0,j0) and ending with (i1,j1) (matrix(i0:i1,j0:j1))
 */
template<typename T>
matrix<T> matrix<T>::getBlock(int i0, int i1, int j0, int j1) const {
	int r = i1 - i0 + 1;
	int c = j1 - j0 + 1;

	if((r < 0) || (c < 0) || (i0 < 0) || (i1 >= _rowCount) || (j0 < 0) || (j1 >= _columnCount))
		throw 0;

	matrix<T> block(r,c);
	for (int i = 0; i < r; ++i)
		memcpy(block.data[i],_data[i+i0]+j0,c*sizeof(T));
	return block;
}

template<typename T>
vector<T> matrix<T>::getColumn(int c, int from, int to) const{
	int r = to-from+1;
	vector<T> ret(r);
	for(int i = from; i < to+1; ++i)
		ret(i) = _data[i][c];
	return ret;
}

template<typename T>
vector<T> matrix<T>::getRow(int r, int from, int to) const{
	int c = to-from+1;
	vector<T> ret(c);
	memcpy(ret._data,_data[r]+from,c*sizeof(T));
	return ret;
}

template<typename T>
matrix<T> decomposeLU(matrix<T> const & m){
	matrix<T> lu(m);
	int n = m.rowCount();
	for(int k = 0; k < n;++k){
		vector<T> col = lu.getColumn(k,0,n-1);
		for(int i = k+1;i<n;++i)
			lu(k,i)/=col(k);
		vector<T> row= lu.getRow(k,0,n-1);
		for(int i = k+1;i < n; ++i)
			for(int j =k+1;j < n; ++j){
				lu(i,j) -=col(i)*row(j);
			}
		//lu(k,k)/=col(k);
	}
	return lu;
}

}

template<typename T>
typename pxx::matrix<T> operator*(T d, typename pxx::matrix<T> const & m) {
	return m.muls(d);
}




#endif /* BLOCK_H_ */
